﻿using System;
using System.Collections.Generic;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
using System.Drawing.Imaging;

namespace Mclib
{
	public enum BlurKernelType
	{
		Exponential,
		Gaussian
	}
	/// <summary>
	/// A class that represents a collection of DataSurface objects.
	/// </summary>
	public class DataSurfaces
	{
		/// <summary>
		/// The DataSurface instances.
		/// </summary>
		public DataSurface[] Surfaces;
		/// <summary>
		/// This is the data defined as a function of surface coordinates (y,x) as Z --> f(y,x)
		/// Z is tabulated at regular intervals over y and x.  Z[row,col] --> Z[y,x]
		/// </summary>
		public double[][,] Z;
		/// <summary>
		/// The value of x corresponding to the first column of Z.
		/// </summary>
		public double Xmin;
		/// <summary>
		/// The value of x corresponding to the last column of Z.
		/// </summary>
		public double Xmax;
		/// <summary>
		/// The value of y corresponding to the first row of Z.
		/// </summary>
		public double Ymin;
		/// <summary>
		/// The value of y corresponding to the last row of Z.
		/// </summary>
		public double Ymax;
		/// <summary>
		/// The blurring level of Zblur.  This is the standard deviation of the Gaussian blur kernel as a proportion of
		/// Z[iSurf].GetLength(0) and Z[iSurf].GetLength(1).
		/// </summary>
		public double BlurLevel
		{
			get
			{
				return blurLevel;
			}
			set
			{
				if (value!=blurLevel)
				{
					blurLevel = value;
					//Surfaces[0].BlurLevel = value;
					for(int iSurf=0; iSurf<Surfaces.Length; iSurf++)
						Surfaces[iSurf].BlurLevel = value;
						//Surfaces[iSurf].SetKern(this.bkType,value,Surfaces[0].blurKern);
				}
			}
		}
		private double blurLevel = 0.0;
		/// <summary>
		/// The shape of blurring kernel.
		/// </summary>
		public BlurKernelType BkType
		{
			get{ return this.bkType; }
			set
			{
				if(value!=this.bkType)
				{
					this.bkType = value;
					Surfaces[0].BkType = value;
					for(int iSurf=1; iSurf<Surfaces.Length; iSurf++)
						Surfaces[iSurf].SetKern(value,blurLevel,Surfaces[0].blurKern);
					if(this.bmp!=null)
						calcBitmap(true);
				}
			}
		}
		private BlurKernelType bkType = BlurKernelType.Exponential;
		/// <summary>
		/// Create a DataSurfaces instance from the specified parameters.
		/// </summary>
		/// <param name="Xmin">The value of x corresponding to the first column of Z.</param>
		/// <param name="Xmax">The value of x corresponding to the last column of Z.</param>
		/// <param name="Ymin">The value of y corresponding to the first row of Z.</param>
		/// <param name="Ymax">The value of y corresponding to the last row of Z.</param>
		/// <param name="Z">This is the data defined as a function of surface coordinates (y,x) as Z --> f(y,x)
		/// Z is tabulated at regular intervals over y and x.  Z[iSurface][row,col] --> Z[iSurface][y,x]</param>
		public DataSurfaces(double Xmin, double Xmax, double Ymin, double Ymax, double[][,] Z)
		{
			this.Xmin = Xmin;
			this.Xmax = Xmax;
			this.Ymin = Ymin;
			this.Ymax = Ymax;
			this.Z = Z;
			Surfaces = new DataSurface[Z.Length];
			double[][] colors = ColorMap.Discrete(Z.Length);
			for(int iSurf = 0; iSurf<Z.Length; iSurf++)
			{
				Surfaces[iSurf] = new DataSurface(Xmin,Xmax,Ymin,Ymax,Z[iSurf]);
				Surfaces[iSurf].Color = colors[iSurf];
			}
		}
		/// <summary>
		/// Create a DataSurface instance from the specified 2D histogram.
		/// </summary>
		/// <param name="h">The 2D histograms.  These must all be defined on the same bins.</param>
		/// <param name="includeTails">Includes the tails of the histogram (outside the range of the axis).</param>
		public DataSurfaces(Histogram2[] h, bool includeTails)
		{
			Surfaces = new DataSurface[h.Length];
			double[][] colors = ColorMap.Discrete(h.Length);
			Z = new double[h.Length][,];
			for(int iSurf=0; iSurf<h.Length; iSurf++)
			{
				Surfaces[iSurf] = new DataSurface(h[iSurf],includeTails);
				Z[iSurf] = Surfaces[iSurf].Z;
				Surfaces[iSurf].Color = colors[iSurf];
			}
			this.Xmin = Surfaces[0].Xmin;
			this.Xmax = Surfaces[0].Xmax;
			this.Ymin = Surfaces[0].Ymin;
			this.Ymax = Surfaces[0].Ymax;
			this.blurLevel = Surfaces[0].BlurLevel;
		}
		/// <summary>
		/// The gamma correction factor used for generating this.Bmp.  Usually a number larger than 1.0 (suggest 2.0).
		/// </summary>
		public double Gamma
		{
			get
			{
				return gamma;
			}
			set
			{
				gamma = value;
				calcBitmap(false);
			}
		}
		private double gamma = 2.0;
		/// <summary>
		/// The bitmap representation of the data surfaces.
		/// </summary>
		public Bitmap Bmp
		{
			get
			{
				if(bmp==null)
					calcBitmap(true);
				return bmp;
			}
		}
		private Bitmap bmp = null;
		private float[,,] bmpCache = null;
		/// <summary>
		/// Pre-calculate the bitmap this.Bmp (before this.Bmp is evaluated for the first time).
		/// </summary>
		public void CalcBitmap()
		{
			calcBitmap(true);
		}
		/// <summary>
		/// Calculates bitmap from the data series.  Used by this.Bmp.
		/// </summary>
		/// <param name="calcCache">If true, this.bmpCache is recalculated (the source of the bitmap data).</param>
		private void calcBitmap(bool calcCache)
		{
			float fr, fg, fb;
			int i,j;
			int iLen = Z[0].GetLength(0);
			int jLen = Z[0].GetLength(1);
			if(bmpCache==null || bmpCache.GetLength(0)<iLen || bmpCache.GetLength(1)<jLen)
			{
				bmpCache = new float[iLen,jLen,3];
				calcCache = true;
			}
			if(calcCache)
			{
				float invRgbMaxFloat;
				double[] rgb = new double[3];
				double v;
				double rgbMax = 0.0;
				int iSurf;
				double[][,] surf = new double[Surfaces.Length][,];
				//	For each surface
				for(iSurf=0; iSurf<Surfaces.Length; iSurf++)
				{
					if(Surfaces[iSurf].Zblur!=null)
						surf[iSurf] = Surfaces[iSurf].Zblur;
					else
						surf[iSurf] = Surfaces[iSurf].Z;
				}
				for(i=0; i<iLen; i++)
					for(j=0; j<jLen; j++)
					{
						//	Zero it out
						rgb[0]=0.0; rgb[1]=0.0; rgb[2]=0.0;
						//	For each surface
						for(iSurf=0; iSurf<Surfaces.Length; iSurf++)
						{
							v = surf[iSurf][i,j];
							rgb[0] += v*Surfaces[iSurf].Color[0];
							rgb[1] += v*Surfaces[iSurf].Color[1];
							rgb[2] += v*Surfaces[iSurf].Color[2];
						}
						//	Running calculation of the maximum value
						rgbMax = Math.Max(rgbMax,Math.Max(rgb[0]+rgb[1],Math.Max(rgb[1]+rgb[2],rgb[2]+rgb[0])));
						//	Store in cache
						bmpCache[i,j,0] = (float)rgb[0];
						bmpCache[i,j,1] = (float)rgb[1];
						bmpCache[i,j,2] = (float)rgb[2];
					}
				//	Normalize cache by maximum value (so that max is 1.0).  Multiplication is faster than division.
				invRgbMaxFloat = 1.0f/(float)rgbMax;
				for(i=0; i<iLen; i++)
					for(j=0; j<jLen; j++)
					{
						fr = bmpCache[i,j,0] * invRgbMaxFloat;
						fg = bmpCache[i,j,1] * invRgbMaxFloat;
						fb = bmpCache[i,j,2] * invRgbMaxFloat;

						bmpCache[i,j,0] = 1.0f-fg-fb;
						bmpCache[i,j,1] = 1.0f-fr-fb;
						bmpCache[i,j,2] = 1.0f-fr-fg;
						/*
						bmpCache[i,j,0] *= invRgbMaxFloat;
						bmpCache[i,j,1] *= invRgbMaxFloat;
						bmpCache[i,j,2] *= invRgbMaxFloat;
						*/
					}
			}
			//	Create bitmap
			if(bmp==null || bmp.Width!=jLen || bmp.Height!=iLen)
				bmp = new Bitmap(jLen,iLen, PixelFormat.Format32bppArgb);
			//	Write pixel information using a fast (unsafe) method.  Safe method is very slow.
			unsafe
			{
				double invGamma = 1.0/gamma;
				BitmapData bd = bmp.LockBits(new Rectangle(0,0,bmp.Width,bmp.Height),ImageLockMode.WriteOnly,PixelFormat.Format32bppArgb);
                byte* ppix;
                int stride = bd.Stride,str;
                int scan0 = bd.Scan0.ToInt32();
				int r,g,b,a;
				float fa,ffa;
				for(i=0; i<iLen; i++)
				{
					str = (iLen-i-1)*stride; // Bitmap is painted upside-down because y-axis is upside-down relative to image axis.
					for(j=0; j<jLen; j++)
					{
						fr = (float)Math.Pow(bmpCache[i,j,0],gamma);
						fg = (float)Math.Pow(bmpCache[i,j,1],gamma);
						fb = (float)Math.Pow(bmpCache[i,j,2],gamma);
						fa = Math.Min(fr,Math.Min(fg,fb));
						ffa = 1.0f/Math.Max(1.0f-fa,0.001f);
						fr = (fr-fa)*ffa;
						fg = (fg-fa)*ffa;
						fb = (fb-fa)*ffa;
						a = 255-(int)(255.9999*fa);
						r = (int)(255.9999*fr);
						g = (int)(255.9999*fg);
						b = (int)(255.9999*fb);
						ppix = (byte*)(str+4*j+scan0);
						* ppix		= (byte)b; // Alpha
						*(ppix+1)	= (byte)g; // Green
						*(ppix+2)	= (byte)r; // Red
						*(ppix+3)	= (byte)a; // Alpha
					}
				}
                bmp.UnlockBits(bd);
			}
		}
	}
	/// <summary>
	/// A class that represents a standard data surface where Z (the value everywhere in the surface)
	/// is a function of x (horizontal coordinate) and y (vertical coordinate).  This class can be used
	/// to draw contour and temperature plots.  This class assumes that sampling of Z is tabulated at
	/// regular intervals over the coordinates x and y.
	/// </summary>
	public class DataSurface
	{
		/// <summary>
		/// This is the data defined as a function of surface coordinates (y,x) as Z --> f(y,x)
		/// Z is tabulated at regular intervals over y and x.  Z[row,col] --> Z[y,x].
		/// The maximum entry is 1.0.
		/// </summary>
		public double[,] Z;
		/// <summary>
		/// The value of x corresponding to the first column of Z.
		/// </summary>
		public double Xmin;
		/// <summary>
		/// The value of x corresponding to the last column of Z.
		/// </summary>
		public double Xmax;
		/// <summary>
		/// The value of y corresponding to the first row of Z.
		/// </summary>
		public double Ymin;
		/// <summary>
		/// The value of y corresponding to the last row of Z.
		/// </summary>
		public double Ymax;
		/// <summary>
		/// The blurring level of Zblur.  This is the standard deviation of the Gaussian blur kernel as a proportion of
		/// Z.GetLength(0) and Z.GetLength(1).  Setting this triggers automatic recalculation of the blur kernel and blurring.
		/// If a bitmap was already created, setting this automaticaly re-paints it.
		/// </summary>
		public double BlurLevel
		{
			get
			{
				return blurLevel;
			}
			set
			{
				if (value!=blurLevel)
				{
					blurLevel = value;
					if (blurLevel <= 0.0)
						Zblur = null;
					else
						this.doBlurring();
					if(this.bmp!=null)
						CalcBitmap();
				}
			}
		}
		private double blurLevel = 0.0;
		/// <summary>
		/// The shape of blurring kernel.  Setting this triggers automatic recalculation of the blur kernel and blurring.
		/// If a bitmap was already created, setting this automaticaly re-paints it.
		/// </summary>
		public BlurKernelType BkType
		{
			get{ return this.bkType; }
			set
			{
				if(value!=this.bkType)
				{
					this.bkType = value;
					this.doBlurring();
					if(this.bmp!=null)
						CalcBitmap();
				}
			}
		}
		private BlurKernelType bkType = BlurKernelType.Gaussian;
		/// <summary>
		/// The blurred surface.  The maximum entry is 1.0.
		/// </summary>
		public double[,] Zblur;
		/// <summary>
		/// The mass of Zblur (the sum of all the entries in Zblur).  Automatically computed when this.BlurLevel is changed (via this.doBlurring()).
		/// </summary>
		public double Sum_Zblur;
		/// <summary>
		/// The blurring kernel computed inside this.doBlurring().
		/// </summary>
		public double[,] blurKern;
		/// <summary>
		/// Create a DataSurface instance from the specified parameters.
		/// </summary>
		/// <param name="Xmin">The value of x corresponding to the first column of Z.</param>
		/// <param name="Xmax">The value of x corresponding to the last column of Z.</param>
		/// <param name="Ymin">The value of y corresponding to the first row of Z.</param>
		/// <param name="Ymax">The value of y corresponding to the last row of Z.</param>
		/// <param name="Z">This is the data defined as a function of surface coordinates (y,x) as Z --> f(y,x)
		/// Z is tabulated at regular intervals over y and x.  Z[row,col] --> Z[y,x]</param>
		public DataSurface(double Xmin, double Xmax, double Ymin, double Ymax, double[,] Z)
		{
			this.Xmin = Xmin;
			this.Xmax = Xmax;
			this.Ymin = Ymin;
			this.Ymax = Ymax;
			this.Z = Z;
		}
		/// <summary>
		/// Create a DataSurface instance from the specified 2D histogram.
		/// </summary>
		/// <param name="h">The 2D histogram.</param>
		/// <param name="includeTails">Includes the tails of the histogram (outside the range of the axis).</param>
		public DataSurface(Histogram2 h, bool includeTails)
		{
			if(includeTails)
			{
				this.Ymin = h.Centers[0][0];
				this.Ymax = h.Centers[0][h.Centers[0].Length-1];
				this.Xmin = h.Centers[1][0];
				this.Xmax = h.Centers[1][h.Centers[1].Length-1];
			}
			else
			{
				this.Ymin = h.LowBinMin[0];
				this.Ymax = h.HighBinMax[0];
				this.Xmin = h.LowBinMin[1];
				this.Xmax = h.HighBinMax[1];
			}
			int max; int total;
			h.Summarize(includeTails, out total, out max);
			this.Z = h.Density(includeTails, max);
			//	Set the blur level and automatically compute stuff.
			this.BlurLevel = 1.0/Math.Sqrt((double)(144+total));
		}
		/// <summary>
		/// Makes a blur kernel and Zblur based on this.BlurLevel.  This is automatically called when this.BlurLevel is set to a new value.
		/// The maximum value of Zblur is automatically set to 1.0.
		/// </summary>
		private void doBlurring()
		{
			//	Dimensions of DataSurface
			int iLen = Z.GetLength(0);
			int jLen = Z.GetLength(1);
			//	Malloc Zblur
			if (Zblur==null || Zblur.GetLength(0)!=iLen || Zblur.GetLength(1)!=jLen)
				Zblur = new double[iLen,jLen];

			//	Indexers of Z, Zblur, and blurKern
			int iz,jz,ik,jk,ii,jj,ikciz,jkcjz;
			//	Convolution bounds on Z spanning the kernel
			int iiMin,iiMax,jjMin,jjMax;
			//	The distance along 'x'
			double ikdsq, jkdsq;
			//	The kernel density at Zblur(x,y) summed across all Z.
			double dens,maxDens=0.0;
			//	Initialize total density to 0.0
			this.Sum_Zblur = 0.0;

			// --------------------------------------------------------------------------
			//	Pre-compute the Gaussian kernel
			// --------------------------------------------------------------------------
			//	Dimensions of blur kernel (go out 3.0 standard deviations on the Gaussian and 4.0 on the exponential)
			double kSpan = 3.0;
			if(bkType==BlurKernelType.Exponential)
				kSpan = 4.0;
			int ikLen = (int)((double)iLen*2.0*kSpan*blurLevel); 
			int jkLen = (int)((double)jLen*2.0*kSpan*blurLevel);
			//	Limit the length if the Gaussian is too big.
			if (ikLen > iLen*2)
				ikLen = iLen*2;
			if (jkLen > jLen*2)
				jkLen = jLen*2;
			//	Make sure the length of the kernel is odd (so that a pixel at the center of the kernel is exactly at the center of the Gaussian).
			if (ikLen%2==0)
				ikLen = ikLen-1;
			if (jkLen%2==0)
				jkLen = jkLen-1;
			//	Get the center coordinates.
			int ikc = ikLen/2; // Center 'i' coordinate
			int jkc = jkLen/2; // Center 'j' coordinate
			//	Precompute a scalar for the kernel
			double si = 1.0/(blurLevel*iLen); // standard deviation is a half scale
			double sj = 1.0/(blurLevel*jLen);
			if(bkType==BlurKernelType.Gaussian)
			{
				si/=Math.Sqrt(2);
				sj/=Math.Sqrt(2);
			}
			double limRadsq = (double)(ikc*ikc*si*si+jkc*jkc*sj*sj)/2.0;
			double radsq;
			//	Malloc kernel
			if (blurKern==null || blurKern.GetLength(0)<ikLen || blurKern.GetLength(1)<jkLen)
				blurKern = new double[ikLen,jkLen];
			//	Compute kernel
			for (ik=0; ik<ikLen; ik++)
			{
				ikdsq = (double)(ik-ikc)*si; // The distance along x
				ikdsq *= ikdsq;
				for (jk=0; jk<jkLen; jk++)
				{
					jkdsq = (double)(jk-jkc)*sj;
					jkdsq *= jkdsq;
					radsq = ikdsq+jkdsq;
					if(radsq>limRadsq)
						blurKern[ik,jk] = 0.0;
					else
					{
						// The distance inside a Gaussian function (the function is scaled to a max of 1).
						if(bkType==BlurKernelType.Gaussian)
							blurKern[ik,jk] = Util.FastExp.Eval(-(radsq));
						else
							blurKern[ik,jk] = Util.FastExp.Eval(-Math.Sqrt(radsq));
					}
				}
			}
			// --------------------------------------------------------------------------

			// --------------------------------------------------------------------------
			//	Convolve
			// --------------------------------------------------------------------------
			//	For each coordinate of Zblur
			for (iz=0; iz<iLen; iz++)
			{
				//	The center of the kernel in Z-coordinates
				ikciz = ikc-iz;
				//	Minimum i-index of Z
				iiMin = iz-ikc;
				if (iiMin<0)
					iiMin = 0;
				//	Maximum i-index of Z
				iiMax = iz+ikc;
				if (iiMax>=iLen)
					iiMax = iLen-1;
				for (jz=0; jz<jLen; jz++)
				{
					//	The center of the kernel in Z-coordinates
					jkcjz = jkc-jz;
					//	Initialize convolution sum to 0.0.
					dens = 0.0;
					//	Minimum j-index of Z
					jjMin = jz-jkc;
					if (jjMin<0)
						jjMin = 0;
					//	Maximum j-index of Z
					jjMax = jz+jkc;
					if (jjMax>=jLen)
						jjMax = jLen-1;
					//	For each i-indez of Z
					for (ii=iiMin; ii<=iiMax; ii++)
					{
						//	Kernel coordinate
						ik = ii+ikciz;
						//	For each j-index of Z
						for (jj=jjMin; jj<=jjMax; jj++)
						{
							//	Increment convolution sum
							dens += blurKern[ik,jj+jkcjz]*Z[ii,jj];
						}
					}
					Zblur[iz,jz] = dens;
					this.Sum_Zblur += dens;
					if (dens >= maxDens)
						maxDens = dens;
				}
			}
			// --------------------------------------------------------------------------

			// --------------------------------------------------------------------------
			//	Normalize to maximum bin size of 1.0
			// --------------------------------------------------------------------------
			si = 1.0/maxDens;
			this.Sum_Zblur*=si;
			for (iz=0; iz<iLen; iz++)
				for (jz=0; jz<jLen; jz++)
					Zblur[iz,jz]*=si;
			// --------------------------------------------------------------------------
		}
		/// <summary>
		/// Allows a pre-computed kernel and the associated parameters to be set without this method
		/// automatically calculating the kernel (via a time-consuming process) which automatically occurs
		/// when this.BkType and this.BlurLevel are set.  After the kernel is set, this method executes
		/// the blurring (convolution) operation.
		/// </summary>
		/// <param name="bkType">The blurring kernel type.</param>
		/// <param name="blurLevel">The blur level used to generate blurKern.</param>
		/// <param name="blurKern">The pre-computed blurring kernel.</param>
		public void SetKern(BlurKernelType bkType, double blurLevel, double[,] blurKern)
		{
			this.blurLevel = blurLevel;
			this.bkType = bkType;
			this.blurKern = blurKern;
		}
		/// <summary>
		/// The color of each surface as [r,g,b] in the range [0,1]
		/// </summary>
		public double[] Color
		{
			get { return color; }
			set
			{
				color = value;
				CalcBitmap();
			}
		}
		private double[] color = new double[] {1.0,0.0,0.0};
		/// <summary>
		/// The gamma correction factor used for generating this.Bmp.  Usually a number larger than 1.0 (suggest 2.0).
		/// </summary>
		public double Gamma
		{
			get
			{
				return gamma;
			}
			set
			{
				gamma = value;
				CalcBitmap();
			}
		}
		private double gamma = 2.0;
		/// <summary>
		/// The bitmap representation of the data surface.
		/// </summary>
		public Bitmap Bmp
		{
			get
			{
				if(bmp==null)
					CalcBitmap();
				return bmp;
			}
		}
		private Bitmap bmp = null;
		/// <summary>
		/// Pre-calculate the bitmap this.Bmp (before this.Bmp is evaluated for the first time).
		/// </summary>
		public void CalcBitmap()
		{
			int i,j;
			double[,] z = Z;
			if(blurLevel>0.0 && Zblur!=null)
				z = Zblur;
			int iLen = Z.GetLength(0);
			int jLen = Z.GetLength(1);
			//	Create bitmap
			if(bmp==null || bmp.Width!=jLen || bmp.Height!=iLen)
				bmp = new Bitmap(jLen,iLen, PixelFormat.Format32bppArgb);
			//	Write pixel information using a fast (unsafe) method.  Safe method is very slow.
			unsafe
			{
				BitmapData bd = bmp.LockBits(new Rectangle(0,0,bmp.Width,bmp.Height),ImageLockMode.WriteOnly,PixelFormat.Format32bppArgb);
                byte* ppix;
                int stride = bd.Stride,str;
                int scan0 = bd.Scan0.ToInt32();
				int r,g,b,a;
				float fr, fg, fb ,fa, ffa;
				double icr = 1.0-Color[0];
				double icg = 1.0-Color[1];
				double icb = 1.0-Color[2];
				double zval;
				for(i=0; i<iLen; i++)
				{
					str = (iLen-i-1)*stride; // Bitmap is painted upside-down because y-axis is upside-down relative to image axis.
					for(j=0; j<jLen; j++)
					{
						zval = z[i,j];
						fr = (float)Math.Pow(1.0-zval*icr,gamma);
						fg = (float)Math.Pow(1.0-zval*icg,gamma);
						fb = (float)Math.Pow(1.0-zval*icb,gamma);
						fa = Math.Min(fr,Math.Min(fg,fb));
						ffa = 1.0f/Math.Max(1.0f-fa,0.001f);
						fr = (fr-fa)*ffa;
						fg = (fg-fa)*ffa;
						fb = (fb-fa)*ffa;
						a = 255-(int)(255.9999*fa);
						r = (int)(255.9999*fr);
						g = (int)(255.9999*fg);
						b = (int)(255.9999*fb);
						ppix = (byte*)(str+4*j+scan0);
						* ppix		= (byte)b; // Alpha
						*(ppix+1)	= (byte)g; // Green
						*(ppix+2)	= (byte)r; // Red
						*(ppix+3)	= (byte)a; // Alpha
					}
				}
                bmp.UnlockBits(bd);
			}
		}
	}
}